import numpy as np
from scipy import linalg, optimize, interpolate
from typing import Dict, List, Tuple, Any
import json
import sys
import traceback

class MathModelingMCP:
    """
    数学建模领域的MCP服务，包含以下领域：
    1. 数值计算算法
    2. 微分方程求解
    3. 数值积分
    4. 优化算法
    5. 插值与拟合算法
    """

    def __init__(self):
        self.methods = {
            "linear_algebra": {
                "gaussian_elimination": self.gaussian_elimination,
                "lu_decomposition": self.lu_decomposition,
                "qr_decomposition": self.qr_decomposition,
                "power_method": self.power_method,
                "jacobi_eigenvalues": self.jacobi_eigenvalues
            },
            "differential_equations": {
                "euler_method": self.euler_method,
                "runge_kutta": self.runge_kutta,
                "finite_difference_pde": self.finite_difference_pde,
                "finite_element": self.finite_element
            },
            "numerical_integration": {
                "trapezoidal_rule": self.trapezoidal_rule,
                "simpsons_rule": self.simpsons_rule,
                "gaussian_quadrature": self.gaussian_quadrature
            },
            "optimization": {
                "simplex_method": self.simplex_method,
                "gradient_descent": self.gradient_descent,
                "newton_optimization": self.newton_optimization,
                "genetic_algorithm": self.genetic_algorithm
            },
            "interpolation_fitting": {
                "lagrange_interpolation": self.lagrange_interpolation,
                "newton_interpolation": self.newton_interpolation,
                "least_squares": self.least_squares,
                "nonlinear_regression": self.nonlinear_regression
            }
        }

    def execute(self, method_path: str, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        执行指定的数学方法
        
        参数:
            method_path: 方法路径，格式为"domain.method_name"
            params: 包含方法所需参数的字典
            
        返回:
            包含执行结果或错误信息的字典
        """
        try:
            domain, method_name = method_path.split('.')
            if domain not in self.methods:
                return {"error": f"Domain '{domain}' not found"}
            
            if method_name not in self.methods[domain]:
                return {"error": f"Method '{method_name}' not found in domain '{domain}'"}
            
            method = self.methods[domain][method_name]
            result = method(params)
            return {"result": result}
            
        except Exception as e:
            return {
                "error": str(e),
                "traceback": traceback.format_exc()
            }

    # 线性代数算法
    def gaussian_elimination(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        高斯消元法求解线性方程组Ax = b
        
        参数:
            A: 系数矩阵 (二维数组)
            b: 常数项向量 (一维数组)
            
        返回:
            x: 方程组的解 (一维数组)
        """
        A = np.array(params['A'], dtype=float)
        b = np.array(params['b'], dtype=float)
        
        if A.shape[0] != A.shape[1]:
            raise ValueError("Coefficient matrix must be square")
        
        if A.shape[0] != b.shape[0]:
            raise ValueError("Matrix and vector dimensions do not match")
        
        n = A.shape[0]
        M = np.hstack([A, b.reshape(-1, 1)])
        for i in range(n):
            # 主元选择
            max_row = i + np.argmax(np.abs(M[i:, i]))
            M[[i, max_row]] = M[[max_row, i]]
            
            # 消元
            for j in range(i+1, n):
                factor = M[j, i] / M[i, i]
                M[j, i:] -= factor * M[i, i:]
        
        # 回代
        x = np.zeros(n)
        for i in range(n-1, -1, -1):
            x[i] = (M[i, -1] - np.sum(M[i, i+1:n] * x[i+1:n])) / M[i, i]
        
        return {
            "solution": x.tolist(),
            "matrix_rank": np.linalg.matrix_rank(A).item(),
            "determinant": np.linalg.det(A).item() if n == m else None
        }

    def lu_decomposition(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        LU分解求解线性方程组Ax = b
        
        参数:
            A: 系数矩阵 (二维数组)
            b: 常数项向量 (一维数组)
            
        返回:
            x: 方程组的解 (一维数组)
            L: 下三角矩阵
            U: 上三角矩阵
        """
        A = np.array(params['A'], dtype=float)
        b = np.array(params['b'], dtype=float)
        
        if A.shape[0] != A.shape[1]:
            raise ValueError("Coefficient matrix must be square")
        
        if A.shape[0] != b.shape[0]:
            raise ValueError("Matrix and vector dimensions do not match")
        
        n = A.shape[0]
        L = np.eye(n)
        U = A.copy()
        
        for i in range(n-1):
            for j in range(i+1, n):
                factor = U[j, i] / U[i, i]
                L[j, i] = factor
                U[j, i:] -= factor * U[i, i:]
        
        # 前向替换
        y = np.zeros(n)
        for i in range(n):
            y[i] = b[i] - np.sum(L[i, :i] * y[:i])
        
        # 回代
        x = np.zeros(n)
        for i in range(n-1, -1, -1):
            x[i] = (y[i] - np.sum(U[i, i+1:n] * x[i+1:n])) / U[i, i]
        
        return {
            "solution": x.tolist(),
            "L": L.tolist(),
            "U": U.tolist(),
            "condition_number": np.linalg.cond(A).item()
        }

    def qr_decomposition(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        QR分解求解线性方程组Ax = b
        
        参数:
            A: 系数矩阵 (二维数组)
            b: 常数项向量 (一维数组)
            
        返回:
            x: 方程组的解 (一维数组)
            Q: 正交矩阵
            R: 上三角矩阵
        """
        A = np.array(params['A'], dtype=float)
        b = np.array(params['b'], dtype=float)
        
        if A.shape[0] != A.shape[1]:
            raise ValueError("Coefficient matrix must be square")
        
        if A.shape[0] != b.shape[0]:
            raise ValueError("Matrix and vector dimensions do not match")
        
        n = A.shape[0]
        Q, R = np.linalg.qr(A)
        
        # 计算Q^T b
        Qtb = np.dot(Q.T, b)
        
        # 回代求解Rx = Qtb
        x = np.zeros(n)
        for i in range(n-1, -1, -1):
            x[i] = (Qtb[i] - np.sum(R[i, i+1:n] * x[i+1:n])) / R[i, i]
        
        return {
            "solution": x.tolist(),
            "Q": Q.tolist(),
            "R": R.tolist(),
            "rank": np.linalg.matrix_rank(A).item()
        }

    def power_method(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        幂法计算矩阵的主特征值和对应特征向量
        
        参数:
            A: 输入矩阵 (二维数组)
            max_iter: 最大迭代次数 (整数, 默认1000)
            tol: 收敛容差 (浮点数, 默认1e-6)
            
        返回:
            eigenvalue: 主特征值
            eigenvector: 对应的特征向量
        """
        A = np.array(params['A'], dtype=float)
        max_iter = params.get('max_iter', 1000)
        tol = params.get('tol', 1e-6)
        
        if A.shape[0] != A.shape[1]:
            raise ValueError("Input matrix must be square")
        
        n = A.shape[0]
        x = np.random.rand(n)
        x = x / np.linalg.norm(x)
        
        for _ in range(max_iter):
            x_new = np.dot(A, x)
            x_new = x_new / np.linalg.norm(x_new)
            
            if np.abs(np.dot(x, x_new) - 1) < tol:
                break
            
            x = x_new
        
        eigenvalue = np.dot(x, np.dot(A, x)) / np.dot(x, x)
        
        return {
            "eigenvalue": eigenvalue.item(),
            "eigenvector": x.tolist(),
            "converged": bool(np.abs(np.dot(x, np.dot(A, x)) - eigenvalue * np.dot(x, x)) < tol)
        }

    def jacobi_eigenvalues(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        雅可比方法计算对称矩阵的所有特征值和特征向量
        
        参数:
            A: 输入对称矩阵 (二维数组)
            max_iter: 最大迭代次数 (整数, 默认1000)
            tol: 收敛容差 (浮点数, 默认1e-6)
            
        返回:
            eigenvalues: 所有特征值 (按降序排列)
            eigenvectors: 对应的特征向量
        """
        A = np.array(params['A'], dtype=float)
        max_iter = params.get('max_iter', 1000)
        tol = params.get('tol', 1e-6)
        
        if A.shape[0] != A.shape[1]:
            raise ValueError("Input matrix must be square")
        
        if not np.allclose(A, A.T):
            raise ValueError("Input matrix must be symmetric")
        
        n = A.shape[0]
        V = np.eye(n)
        A_diag = np.diag(A).copy()
        
        for _ in range(max_iter):
            off_diag_norm = np.sqrt(np.sum(np.triu(A, 1)**2))
            
            if off_diag_norm < tol:
                break
            
            for p in range(n):
                for q in range(p+1, n):
                    if abs(A[p, q]) > 1e-8:
                        theta = (A[q, q] - A[p, p]) / (2 * A[p, q])
                        t = 1 / (abs(theta) + np.sqrt(1 + theta**2))
                        if theta < 0:
                            t = -t
                        
                        c = 1 / np.sqrt(1 + t**2)
                        s = t * c
                        
                        # 更新矩阵A
                        Apq = A[p, q]
                        App = A[p, p]
                        Aqq = A[q, q]
                        
                        A[p, q] = 0
                        A[q, p] = 0
                        A[p, p] = App - t * Apq
                        A[q, q] = Aqq + t * Apq
                        
                        for r in range(n):
                            if r != p and r != q:
                                Artmp = A[r, p]
                                A[r, p] = Artmp * c - A[r, q] * s
                                A[p, r] = A[r, p]
                                A[r, q] = A[r, q] * c + Artmp * s
                                A[q, r] = A[r, q]
                        
                        # 更新特征向量
                        for r in range(n):
                            Vrtmp = V[r, p]
                            V[r, p] = Vrtmp * c - V[r, q] * s
                            V[r, q] = V[r, q] * c + Vrtmp * s
        
        eigenvalues = np.diag(A)
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx].tolist()
        eigenvectors = V[:, idx].tolist()
        
        return {
            "eigenvalues": eigenvalues,
            "eigenvectors": eigenvectors,
            "converged": bool(np.all(np.abs(np.triu(A, 1)) < tol))
        }

    # 微分方程求解算法
    def euler_method(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        欧拉法求解常微分方程初值问题
        
        参数:
            dy_dx: 微分方程函数，形式为f(x, y)
            x0: 初始x值
            y0: 初始y值
            x_end: 结束x值
            h: 步长
            
        返回:
            x_values: x值列表
            y_values: 对应的y值列表
        """
        dy_dx_str = params['dy_dx']
        dy_dx = eval(f"lambda x, y: {dy_dx_str}")
        x0 = float(params['x0'])
        y0 = float(params['y0'])
        x_end = float(params['x_end'])
        h = float(params['h'])
        
        if h <= 0:
            raise ValueError("Step size must be positive")
        
        n = int((x_end - x0) / h) + 1
        x = np.linspace(x0, x_end, n)
        y = np.zeros(n)
        y[0] = y0
        
        for i in range(n-1):
            y[i+1] = y[i] + h * dy_dx(x[i], y[i])
        
        return {
            "x_values": x.tolist(),
            "y_values": y.tolist(),
            "step_size": h,
            "num_steps": n-1
        }

    def runge_kutta(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        经典四阶龙格-库塔法求解常微分方程初值问题
        
        参数:
            dy_dx: 微分方程函数，形式为f(x, y)
            x0: 初始x值
            y0: 初始y值
            x_end: 结束x值
            h: 步长
            
        返回:
            x_values: x值列表
            y_values: 对应的y值列表
        """
        dy_dx_str = params['dy_dx']
        dy_dx = eval(f"lambda x, y: {dy_dx_str}")
        x0 = float(params['x0'])
        y0 = float(params['y0'])
        x_end = float(params['x_end'])
        h = float(params['h'])
        
        if h <= 0:
            raise ValueError("Step size must be positive")
        
        n = int((x_end - x0) / h) + 1
        x = np.linspace(x0, x_end, n)
        y = np.zeros(n)
        y[0] = y0
        
        for i in range(n-1):
            k1 = h * dy_dx(x[i], y[i])
            k2 = h * dy_dx(x[i] + h/2, y[i] + k1/2)
            k3 = h * dy_dx(x[i] + h/2, y[i] + k2/2)
            k4 = h * dy_dx(x[i] + h, y[i] + k3)
            
            y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6

    def newton_cotes(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        牛顿-柯特斯公式进行数值积分
        
        参数:
            func: 要积分的函数，接受一个参数x并返回数值
            a: 积分下限
            b: 积分上限
            n: 分割区间数
            method: 使用的方法 ('trapezoidal' 梯形法则, 'simpson' 辛普森法则)
            
        返回:
            integral: 积分结果
            error_estimate: 误差估计
        """
        func = params['func']
        a = params['a']
        b = params['b']
        n = params['n']
        method = params.get('method', 'simpson')
        
        h = (b - a) / n
        x = np.linspace(a, b, n+1)
        y = np.array([func(xi) for xi in x])
        
        if method == 'trapezoidal':
            # 梯形公式
            integral = h * (0.5*y[0] + 0.5*y[-1] + np.sum(y[1:-1]))
            # 误差估计（基于余项）
            max_second_derivative = params.get('max_second_derivative', 1e6)  # 默认最大二阶导数
            error_estimate = (b - a)**3 * max_second_derivative / (12 * n**2)
        elif method == 'simpson':
            # 辛普森公式（要求n为偶数）
            if n % 2 != 0:
                raise ValueError("辛普森法则要求分割区间数为偶数")
            
            integral = h/3 * (y[0] + y[-1] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-2:2]))
            # 误差估计（基于余项）
            max_fourth_derivative = params.get('max_fourth_derivative', 1e6)  # 默认最大四阶导数
            error_estimate = (b - a)**5 * max_fourth_derivative / (180 * n**4)
        else:
            raise ValueError(f"未知方法: {method}")
        
        return {"integral": integral, "error_estimate": error_estimate}

    def finite_element(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        有限元法求解偏微分方程
        
        参数:
            pde_type: 偏微分方程类型 ('poisson' 泊松方程, 'helmholtz' 赫姆霍兹方程)
            boundary_conditions: 边界条件，格式为:
                {
                    'left': {'type': 'dirichlet', 'value': 0},
                    'right': {'type': 'dirichlet', 'value': 0},
                    'top': {'type': 'dirichlet', 'value': 0},
                    'bottom': {'type': 'dirichlet', 'value': 0}
                }
            source_term: 源项函数（对于泊松方程）
            wave_number: 波数（对于赫姆霍兹方程）
            x_range: x轴范围 [x_start, x_end]
            y_range: y轴范围 [y_start, y_end]
            num_elements_x: x轴单元数
            num_elements_y: y轴单元数
            element_type: 单元类型 ('linear_triangular' 线性三角形, 'bilinear_rectangular' 双线性矩形)
            
        返回:
            solution: 解的二维数组
            mesh: 网格信息
        """
        pde_type = params['pde_type']
        boundary_conditions = params['boundary_conditions']
        x_range = params.get('x_range', [0, 1])
        y_range = params.get('y_range', [0, 1])
        num_elements_x = params['num_elements_x']
        num_elements_y = params['num_elements_y']
        element_type = params.get('element_type', 'bilinear_rectangular')
        
        # 创建网格
        if element_type == 'linear_triangular':
            # 三角形单元网格
            num_nodes_x = num_elements_x + 1
            num_nodes_y = num_elements_y + 1
            
            # 节点坐标
            x = np.linspace(x_range[0], x_range[1], num_nodes_x)
            y = np.linspace(y_range[0], y_range[1], num_nodes_y)
            X, Y = np.meshgrid(x, y)
            
            # 元素连接
            elements = []
            for i in range(num_elements_x):
                for j in range(num_elements_y):
                    node1 = j*num_nodes_x + i
                    node2 = node1 + 1
                    node3 = node1 + num_nodes_x
                    elements.append([node1, node2, node3])
                    
                    node4 = node3
                    node5 = node2
                    node6 = node3 + 1
                    elements.append([node4, node5, node6])
            
    elif element_type == 'bilinear_rectangular':
            # 双线性矩形单元网格
            num_nodes_x = num_elements_x + 1
            num_nodes_y = num_elements_y + 1
            
            # 节点坐标
            x = np.linspace(x_range[0], x_range[1], num_nodes_x)
            y = np.linspace(y_range[0], y_range[1], num_nodes_y)
            X, Y = np.meshgrid(x, y)
            
            # 元素连接
            elements = []
            for i in range(num_elements_x):
                for j in range(num_elements_y):
                    node1 = j*num_nodes_x + i
                    node2 = node1 + 1
                    node3 = node1 + num_nodes_x + 1
                    node4 = node1 + num_nodes_x
                    elements.append([node1, node2, node3, node4])
            
        # 组装刚度矩阵和载荷向量
        num_nodes = num_nodes_x * num_nodes_y
        K = np.zeros((num_nodes, num_nodes))
        F = np.zeros(num_nodes)
        
        # 计算单元刚度矩阵和载荷向量
        for elem in elements:
            if element_type == 'linear_triangular':
                # 三角形单元计算
                coords = np.array([[X[i], Y[j]] for i, j in [(elem[0]%num_nodes_x, elem[0]//num_nodes_x),         (elem[1]%num_nodes_x, elem[1]//num_nodes_x), (elem[2]%num_nodes_x, elem[2]//num_nodes_x)]])

    def gaussian_quadrature(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        高斯积分法进行数值积分
        
        参数:
            func: 要积分的函数，接受一个参数x并返回数值
            a: 积分下限
            b: 积分上限
            n: 使用的节点数
            
        返回:
            integral: 积分结果
            nodes: 使用的节点
            weights: 对应的权重
            error_estimate: 误差估计
        """
        func = params['func']
        a = params['a']
        b = params['b']
        n = params['n']
        
        # 获取高斯-勒让德节点和权重
        nodes, weights = np.polynomial.legendre.leggauss(n)
        
        # 将积分区间[-1, 1]变换到[a, b]
        transformed_nodes = (b - a)/2 * nodes + (a + b)/2
        transformed_weights = (b - a)/2 * weights
        
        # 计算积分
        integral = np.sum(transformed_weights * np.array([func(x) for x in transformed_nodes]))
        
        # 误差估计（基于余项）
        # 对于高斯积分，误差与2n阶导数有关
        max_2n_derivative = params.get('max_2n_derivative', 1e6)  # 默认最大2n阶导数
        error_estimate = (b - a)**(2*n+1) / (2**2*n * (2*n+1) * (np.math.factorial(n))**2) * max_2n_derivative
        
        return {
            "integral": integral, 
            "nodes": transformed_nodes.tolist(), 
            "weights": transformed_weights.tolist(),
            "error_estimate": error_estimate
        }

    def simplex_method(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        单纯形法求解线性规划问题
        
        参数:
            c: 目标函数系数向量
            A: 约束矩阵
            b: 约束条件右侧值向量
            maximize: 是否是最大化问题（默认为True）
            
        返回:
            solution: 解向量
            objective_value: 目标函数值
            success: 是否成功找到最优解
            iterations: 迭代次数
        """
        c = np.array(params['c'], dtype=float)
        A = np.array(params['A'], dtype=float)
        b = np.array(params['b'], dtype=float)
        maximize = params.get('maximize', True)
        
        if A.shape[0] != b.shape[0]:
            raise ValueError("约束矩阵行数必须与右侧值向量维度相同")
        
        if A.shape[1] != c.shape[0]:
            raise ValueError("约束矩阵列数必须与目标函数系数向量维度相同")
        
        # 处理最大化问题（转换为最小化问题）
        if maximize:
            c = -c
        
        # 添加松弛变量
        m, n = A.shape
        A_augmented = np.hstack([A, np.eye(m)])
        c_augmented = np.hstack([c, np.zeros(m)])
        
        # 初始化单纯形表
        tableau = np.vstack([
            np.hstack([A_augmented, b.reshape(-1, 1)]),
            np.hstack([-c_augmented, 0])
        ])
        
        iterations = 0
        success = False
        
        while True:
            iterations += 1
            
            # 检查是否达到最优条件
            if np.all(tableau[-1, :-1] >= 0):
                success = True
                break
            
            # 选择进入基的变量
            entering_col = np.argmin(tableau[-1, :-1])
            
            # 检查是否有正系数（否则问题无界）
            if np.all(tableau[:-1, entering_col] <= 0):
                success = False
                break
            
            # 选择离开基的变量
            ratios = []
            for i in range(m):
                if tableau[i, entering_col] > 0:
                    ratios.append((tableau[i, -1] / tableau[i, entering_col], i))
            
            leaving_row = min(ratios, key=lambda x: x[0])[1]
            
            # 更新单纯形表
            pivot = tableau[leaving_row, entering_col]
            tableau[leaving_row, :] /= pivot
            
            for i in range(m + 1):
                if i != leaving_row:
                    factor = tableau[i, entering_col]
                    tableau[i, :] -= factor * tableau[leaving_row, :]
        
        # 提取解
        solution = np.zeros(n + m)
        for i in range(m):
            col = np.where(np.abs(tableau[i, :-1]) > 1e-10)[0]
            if len(col) == 1:
                solution[col[0]] = tableau[i, -1]
        
        # 截断非决策变量部分
        solution = solution[:n]
        objective_value = np.dot(c, solution)
        
        return {
            "solution": solution.tolist(),
            "objective_value": objective_value.item() if maximize else -objective_value.item(),
            "success": success,
            "iterations": iterations
        }

    def gradient_descent(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        梯度下降法求解非线性优化问题
        
        参数:
            func: 要优化的目标函数，接受一个参数x并返回数值
            grad_func: 目标函数的梯度函数，接受一个参数x并返回梯度向量
            x0: 初始解向量
            learning_rate: 学习率（步长）
            max_iter: 最大迭代次数
            tol: 收敛容差
            
        返回:
            solution: 最优解向量
            objective_value: 目标函数值
            iterations: 迭代次数
            converged: 是否收敛
        """
        func = params['func']
        grad_func = params['grad_func']
        x0 = np.array(params['x0'], dtype=float)
        learning_rate = params.get('learning_rate', 0.01)
        max_iter = params.get('max_iter', 1000)
        tol = params.get('tol', 1e-6)
        
        x = x0.copy()
        history = []
        
        for iteration in range(max_iter):
            gradient = grad_func(x)
            
            # 更新解
            x_new = x - learning_rate * gradient
            
            # 计算变化量
            change = np.linalg.norm(x_new - x)
            x = x_new
            
            # 记录目标函数值
            objective_value = func(x)
            history.append((iteration, x.copy(), objective_value, change))
            
            # 检查收敛条件
            if change < tol:
                converged = True
                break
            
            if iteration == max_iter - 1:
                converged = False
        
        return {
            "solution": x.tolist(),
            "objective_value": objective_value.item(),
            "iterations": iteration + 1,
            "converged": converged,
            "history": [(iter_, x_.tolist(), obj_val, change_) for iter_, x_, obj_val, change_ in history]
        }

    def newton_method(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        牛顿法求解非线性优化问题
        
        参数:
            func: 要优化的目标函数，接受一个参数x并返回数值
            grad_func: 目标函数的梯度函数，接受一个参数x并返回梯度向量
            hess_func: 目标函数的Hessian矩阵函数，接受一个参数x并返回Hessian矩阵
            x0: 初始解向量
            max_iter: 最大迭代次数
            tol: 收敛容差
            
        返回:
            solution: 最优解向量
            objective_value: 目标函数值
            iterations: 迭代次数
            converged: 是否收敛
            history: 优化过程的历史记录
        """
        func = params['func']
        grad_func = params['grad_func']
        hess_func = params['hess_func']
        x0 = np.array(params['x0'], dtype=float)
        max_iter = params.get('max_iter', 1000)
        tol = params.get('tol', 1e-6)
        
        x = x0.copy()
        history = []
        
        for iteration in range(max_iter):
            gradient = grad_func(x)
            hessian = hess_func(x)
            
            # 计算Hessian矩阵的逆
            try:
                hessian_inv = np.linalg.inv(hessian)
            except np.linalg.LinAlgError:
                raise ValueError("Hessian矩阵不可逆")
            
            # 更新解
            x_new = x - np.dot(hessian_inv, gradient)
            
            # 计算变化量
            change = np.linalg.norm(x_new - x)
            x = x_new
            
            # 记录目标函数值
            objective_value = func(x)
            history.append((iteration, x.copy(), objective_value, change))
            
            # 检查收敛条件
            if change < tol:
                converged = True
                break
            
            if iteration == max_iter - 1:
                converged = False
        
        return {
            "solution": x.tolist(),
            "objective_value": objective_value.item(),
            "iterations": iteration + 1,
            "converged": converged,
            "history": [(iter_, x_.tolist(), obj_val, change_) for iter_, x_, obj_val, change_ in history]
        }

    def genetic_algorithm(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        遗传算法求解非线性优化问题
        
        参数:
            func: 要优化的目标函数，接受一个参数x并返回数值
            bounds: 变量的取值范围，形如[(low1, high1), (low2, high2), ...]
            population_size: 种群大小
            generations: 迭代代数
            crossover_rate: 交叉概率
            mutation_rate: 变异概率
            elitism: 是否使用精英保留策略（默认为True）
            maximize: 是否是最大化问题（默认为False）
            
        返回:
            solution: 最优解向量
            objective_value: 目标函数值
            generation: 达到最优解的代数
            best_fitness_history: 每一代的最佳适应度值
            average_fitness_history: 每一代的平均适应度值
        """
        import random
        import numpy as np
        
        func = params['func']
        bounds = params['bounds']
        population_size = params.get('population_size', 50)
        generations = params.get('generations', 100)
        crossover_rate = params.get('crossover_rate', 0.8)
        mutation_rate = params.get('mutation_rate', 0.1)
        elitism = params.get('elitism', True)
        maximize = params.get('maximize', False)
        
        num_vars = len(bounds)
        
        # 初始化种群
        def initialize_population():
            population = []
            for _ in range(population_size):
                individual = [random.uniform(low, high) for low, high in bounds]
                population.append(individual)
            return population
        
        # 计算适应度
        def calculate_fitness(population):
            fitness_values = []
            for individual in population:
                fitness = func(individual)
                if maximize:
                    fitness_values.append(fitness)
                else:
                    fitness_values.append(-fitness)  # 最小化问题转换为最大化问题
            return np.array(fitness_values)
        
        # 选择操作（轮盘赌选择）
        def selection(population, fitness_values):
            # 归一化适应度
            min_fitness = np.min(fitness_values)
            if min_fitness < 0:
                fitness_values += abs(min_fitness)
            elif min_fitness == 0:
                fitness_values += 1e-10  # 避免除以零
            
            total_fitness = np.sum(fitness_values)
            probabilities = fitness_values / total_fitness
            
            selected_indices = np.random.choice(range(population_size), size=population_size, p=probabilities)
            selected_population = [population[i] for i in selected_indices]
            return selected_population
        
        # 交叉操作（单点交叉）
        def crossover(population):
            new_population = []
            for i in range(0, population_size, 2):
                parent1 = population[i]
                parent2 = population[i+1] if i+1 < population_size else population[0]
                
                if random.random() < crossover_rate:
                    # 单点交叉
                    crossover_point = random.randint(1, num_vars-1)
                    child1 = parent1[:crossover_point] + parent2[crossover_point:]
                    child2 = parent2[:crossover_point] + parent1[crossover_point:]
                else:
                    child1, child2 = parent1, parent2
                
                new_population.extend([child1, child2])
            return new_population
        
        # 变异操作（随机扰动）
        def mutation(population):
            for i in range(len(population)):
                for j in range(num_vars):
                    if random.random() < mutation_rate:
                        # 在变量范围内进行随机扰动
                        perturbation = random.uniform(-0.1, 0.1) * (bounds[j][1] - bounds[j][0])
                        population[i][j] += perturbation
                        # 确保变异后的值仍在边界范围内
                        population[i][j] = max(bounds[j][0], min(bounds[j][1], population[i][j]))
            return population
        
        # 主程序
        population = initialize_population()
        best_fitness_history = []
        average_fitness_history = []
        best_solution = None
        best_fitness = float('-inf') if maximize else float('inf')
        best_generation = 0
        
        for generation in range(generations):
            fitness_values = calculate_fitness(population)
            
            # 记录最佳和平均适应度
            current_best_fitness = np.max(fitness_values) if maximize else -np.min(fitness_values)
            current_average_fitness = np.mean(fitness_values) if maximize else -np.mean(fitness_values)
            best_fitness_history.append(current_best_fitness)
            average_fitness_history.append(current_average_fitness)
            
            # 更新全局最佳解
            if maximize:
                if current_best_fitness > best_fitness:
                    best_fitness = current_best_fitness
                    best_solution = population[np.argmax(fitness_values)]
                    best_generation = generation
            else:
                if current_best_fitness < best_fitness:
                    best_fitness = current_best_fitness
                    best_solution = population[np.argmin(fitness_values)]
                    best_generation = generation
            
            # 精英保留
            if elitism and best_solution is not None:
                new_population = selection(population, fitness_values)
                new_population = crossover(new_population)
                new_population = mutation(new_population)
                
                # 替换最差个体
                worst_index = np.argmin(fitness_values) if maximize else np.argmax(fitness_values)
                new_population[worst_index] = best_solution
                population = new_population
            else:
                population = mutation(crossover(selection(population, fitness_values)))
        
        return {
            "solution": best_solution,
            "objective_value": best_fitness if maximize else -best_fitness,
            "generation": best_generation,
            "best_fitness_history": best_fitness_history,
            "average_fitness_history": average_fitness_history
        }

    def lagrange_interpolation(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        拉格朗日插值法
        
        参数:
            x_data: 已知数据点的x坐标列表
            y_data: 已知数据点的y坐标列表
            x_interp: 要插值的x坐标
            
        返回:
            y_interp: 插值得到的y值
            polynomial: 插值多项式的字符串表示
        """
        import numpy as np
        from sympy import symbols, expand
        
        x_data = np.array(params['x_data'], dtype=float)
        y_data = np.array(params['y_data'], dtype=float)
        x_interp = float(params['x_interp'])
        
        n = len(x_data)
        if n != len(y_data):
            raise ValueError("x_data和y_data的长度必须相同")
        
        # 构造拉格朗日基函数
        x = symbols('x')
        basis_functions = []
        for i in range(n):
            numerator = 1
            denominator = 1
            for j in range(n):
                if i != j:
                    numerator *= (x - x_data[j])
                    denominator *= (x_data[i] - x_data[j])
            basis_functions.append(numerator / denominator)
        
        # 构造插值多项式
        interpolation_poly = 0
        for i in range(n):
            interpolation_poly += y_data[i] * basis_functions[i]
        
        # 计算插值点的值
        y_interp = float(interpolation_poly.subs(x, x_interp))
        
        return {
            "y_interp": y_interp,
            "polynomial": str(expand(interpolation_poly))
        }

    def newton_interpolation(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        牛顿插值法
        
        参数:
            x_data: 已知数据点的x坐标列表
            y_data: 已知数据点的y坐标列表
            x_interp: 要插值的x坐标
            
        返回:
            y_interp: 插值得到的y值
            polynomial: 插值多项式的字符串表示
            divided_differences_table: 差商表
        """
        import numpy as np
        from sympy import symbols, expand
        
        x_data = np.array(params['x_data'], dtype=float)
        y_data = np.array(params['y_data'], dtype=float)
        x_interp = float(params['x_interp'])
        
        n = len(x_data)
        if n != len(y_data):
            raise ValueError("x_data和y_data的长度必须相同")
        
        # 构造差商表
        divided_differences = np.zeros((n, n))
        divided_differences[:,0] = y_data  # 0阶差商就是y值
        
        for j in range(1, n):  # 计算j阶差商
            for i in range(n - j):
                divided_differences[i,j] = (divided_differences[i+1,j-1] - divided_differences[i,j-1]) / (x_data[i+j] - x_data[i])
        
        # 构造插值多项式
        x = symbols('x')
        interpolation_poly = divided_differences[0,0]
        product_term = 1  # 存储(x-x0)(x-x1)...(x-xk)项
        
        for i in range(1, n):
            product_term *= (x - x_data[i-1])
            interpolation_poly += divided_differences[0,i] * product_term
        
        # 计算插值点的值
        y_interp = float(interpolation_poly.subs(x, x_interp))
        
        return {
            "y_interp": y_interp,
            "polynomial": str(expand(interpolation_poly)),
            "divided_differences_table": divided_differences.tolist()
        }

    def linear_least_squares(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        线性最小二乘法拟合数据点
        
        参数:
            x_data: 已知数据点的x坐标列表
            y_data: 已知数据点的y坐标列表
            degree: 拟合多项式的阶数
            
        返回:
            coefficients: 拟合多项式的系数列表（从低阶到高阶）
            residuals: 残差平方和
            equation: 拟合多项式的字符串表示
            predicted_values: 拟合曲线在输入点处的预测值
        """
        import numpy as np
        from sympy import symbols, expand
        
        x_data = np.array(params['x_data'], dtype=float)
        y_data = np.array(params['y_data'], dtype=float)
        degree = params.get('degree', 1)
        
        n = len(x_data)
        if n != len(y_data):
            raise ValueError("x_data和y_data的长度必须相同")
        if n <= degree:
            raise ValueError("数据点数量必须大于多项式阶数")
        
        # 构造设计矩阵
        X = np.vander(x_data, degree + 1, increasing=True)  # 生成范德蒙矩阵
        
        # 使用正规方程求解最小二乘问题
        coefficients = np.linalg.solve(X.T @ X, X.T @ y_data)
        
        # 计算拟合值和残差
        predicted_values = X @ coefficients
        residuals = np.sum((y_data - predicted_values)**2)
        
        # 构造多项式表达式
        x = symbols('x')
        polynomial = 0
        for i in range(degree + 1):
            polynomial += coefficients[i] * x**i
        
        return {
            "coefficients": coefficients.tolist(),
            "residuals": float(residuals),
            "equation": str(expand(polynomial)),
            "predicted_values": predicted_values.tolist()
        }

    def nonlinear_least_squares(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        非线性最小二乘法拟合数据点
        
        参数:
            x_data: 已知数据点的x坐标列表
            y_data: 已知数据点的y坐标列表
            func: 要拟合的非线性函数，形式为 f(x, *params)
            p0: 初始参数猜测值列表
            max_iter: 最大迭代次数
            tol: 收敛容忍度
            
        返回:
            parameters: 拟合得到的参数值列表
            covariance: 参数估计的协方差矩阵
            residuals: 残差平方和
            predicted_values: 拟合曲线在输入点处的预测值
            success: 是否成功收敛
            message: 优化过程的信息
        """
        import numpy as np
        from scipy.optimize import least_squares
        
        x_data = np.array(params['x_data'], dtype=float)
        y_data = np.array(params['y_data'], dtype=float)
        func = params['func']  # 非线性函数
        p0 = np.array(params['p0'], dtype=float)  # 初始参数猜测
        max_iter = params.get('max_iter', 100)
        tol = params.get('tol', 1e-8)
        
        n = len(x_data)
        if n != len(y_data):
            raise ValueError("x_data和y_data的长度必须相同")
        if n <= len(p0):
            raise ValueError("数据点数量必须大于参数数量")
        
        # 定义残差函数
        def residuals_func(p):
            return func(x_data, *p) - y_data
        
        # 使用Levenberg-Marquardt算法
        result = least_squares(residuals_func, p0, method='lm', max_nfev=max_iter, ftol=tol)
        
        # 计算预测值
        predicted_values = func(x_data, *result.x)
        
        return {
            "parameters": result.x.tolist(),
            "covariance": result.jac.tolist() if result.jac is not None else None,
            "residuals": float(np.sum((y_data - predicted_values)**2)),
            "predicted_values": predicted_values.tolist(),
            "success": result.success,
            "message": result.message
        }

    def polynomial_interpolation_and_fitting(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        多项式插值与拟合综合方法
        
        参数:
            x_data: 已知数据点的x坐标列表
            y_data: 已知数据点的y坐标列表
            method: 插值/拟合方法 ('lagrange', 'newton', 'linear_least_squares', 'nonlinear_least_squares')
            degree: 拟合多项式的阶数（仅对线性最小二乘法有效）
            func: 非线性函数（仅对非线性最小二乘法有效）
            p0: 初始参数猜测值列表（仅对非线性最小二乘法有效）
            max_iter: 最大迭代次数（仅对非线性最小二乘法有效）
            tol: 收敛容忍度（仅对非线性最小二乘法有效）
        
        返回:
            根据选择的方法返回相应的结果字典
        """
        import numpy as np
        
        method = params['method']
        x_data = np.array(params['x_data'], dtype=float)
        y_data = np.array(params['y_data'], dtype=float)
        
        if method == 'lagrange':
            # 使用拉格朗日插值
            return self.lagrange_interpolation({
                'x_data': x_data,
                'y_data': y_data,
                'x_interp': x_data  # 返回所有插值点的结果
            })
        
        elif method == 'newton':
            # 使用牛顿插值
            return self.newton_interpolation({
                'x_data': x_data,
                'y_data': y_data,
                'x_interp': x_data  # 返回所有插值点的结果
            })
        
        elif method == 'linear_least_squares':
            # 使用线性最小二乘法拟合
            degree = params.get('degree', 1)
            return self.linear_least_squares({
                'x_data': x_data,
                'y_data': y_data,
                'degree': degree
            })
        
        elif method == 'nonlinear_least_squares':
            # 使用非线性最小二乘法拟合
            if 'func' not in params or 'p0' not in params:
                raise ValueError("对于非线性最小二乘法，必须提供func和p0参数")
            
            func = params['func']
            p0 = np.array(params['p0'], dtype=float)
            max_iter = params.get('max_iter', 100)
            tol = params.get('tol', 1e-8)
            
            return self.nonlinear_least_squares({
                'x_data': x_data,
                'y_data': y_data,
                'func': func,
                'p0': p0,
                'max_iter': max_iter,
                'tol': tol
            })
        
        else:
            raise ValueError(f"不支持的方法类型: {method}")

    def exponential_fit(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        指数函数拟合数据点（形如 y = a * e^(b*x)）
        
        参数:
            x_data: 已知数据点的x坐标列表
            y_data: 已知数据点的y坐标列表
            p0: 初始参数猜测值列表 [a0, b0]
            max_iter: 最大迭代次数
            tol: 收敛容忍度
            
        返回:
            parameters: 拟合得到的参数值列表 [a, b]
            covariance: 参数估计的协方差矩阵
            residuals: 残差平方和
            predicted_values: 拟合曲线在输入点处的预测值
            success: 是否成功收敛
            message: 优化过程的信息
        """
        import numpy as np
        from scipy.optimize import curve_fit
        
        x_data = np.array(params['x_data'], dtype=float)
        y_data = np.array(params['y_data'], dtype=float)
        p0 = np.array(params['p0'], dtype=float) if 'p0' in params else None  # 初始参数猜测
        max_iter = params.get('max_iter', 100)
        tol = params.get('tol', 1e-8)
        
        n = len(x_data)
        if n != len(y_data):
            raise ValueError("x_data和y_data的长度必须相同")
        if p0 is not None and n <= len(p0):
            raise ValueError("数据点数量必须大于参数数量")
        
        # 定义指数函数
        def exp_func(x, a, b):
            return a * np.exp(b * x)
        
        # 使用curve_fit进行非线性拟合
        try:
            popt, pcov = curve_fit(exp_func, x_data, y_data, p0=p0, maxfev=max_iter, ftol=tol)
        except RuntimeError as e:
            return {
                "parameters": None,
                "covariance": None,
                "residuals": None,
                "predicted_values": None,
                "success": False,
                "message": str(e)
            }
        
        # 计算预测值
        predicted_values = exp_func(x_data, *popt)
        
        return {
            "parameters": popt.tolist(),
            "covariance": pcov.tolist() if pcov is not None else None,
            "residuals": float(np.sum((y_data - predicted_values)**2)),
            "predicted_values": predicted_values.tolist(),
            "success": True,
            "message": "Optimization terminated successfully."
        }
    coords = [(elem[1]%num_nodes_x, elem[1]//num_nodes_x),
              (elem[2]%num_nodes_x, elem[2]//num_nodes_x)]
    # 计算面积和形状函数导数
    A = 0.5 * abs((coords[1,0]-coords[0,0])*(coords[2,1]-coords[0,1]) - 
                  (coords[2,0]-coords[0,0])*(coords[1,1]-coords[0,1]))

    B = np.array([
        [coords[1,1]-coords[0,1], coords[2,1]-coords[0,1]],
        [coords[0,0]-coords[1,0], coords[0,0]-coords[2,0]]
    ]) / (2*A)
                
    # 单元刚度矩阵
    Ke = A * np.dot(B.T, B)
                
    # 单元载荷向量（假设源项为常数）
    if 'source_term' in params and pde_type == 'poisson':
        source = params['source_term']
        fe = A * source / 3 * np.ones(3)
    else:
        fe = np.zeros(3)
            
elif element_type == 'bilinear_rectangular':

                # 双线性矩形单元计算
                coords = np.array([[X[i], Y[j]] for i, j in [(elem[0]%num_nodes_x, elem[0]//num_nodes_x),

                                                              (elem[1]%num_nodes_x, elem[1]//num_nodes_x), 
                                                              (elem[2]%num_nodes_x, elem[2]//num_nodes_x), 
                                                              (elem[3]%num_nodes_x, elem[3]//num_nodes_x)]])
                
                # 计算雅可比矩阵和形状函数导数
                jac = np.array([
                    [coords[1,0]-coords[0,0], coords[3,0]-coords[0,0]],
                    [coords[1,1]-coords[0,1], coords[3,1]-coords[0,1]]
                ])
                det_jac = np.linalg.det(jac)
                inv_jac = np.linalg.inv(jac)
                
                # 高斯积分点（2x2积分）
                gauss_points = [
                    (-1/np.sqrt(3), -1/np.sqrt(3)),
                    (1/np.sqrt(3), -1/np.sqrt(3)),
                    (1/np.sqrt(3), 1/np.sqrt(3)),
                    (-1/np.sqrt(3), 1/np.sqrt(3))
                ]
                weights = [1, 1, 1, 1]
                
                # 计算单元刚度矩阵
                Ke = np.zeros((4, 4))
                for gp, w in zip(gauss_points, weights):
                    xi, eta = gp
                    
                    # 形状函数导数
                    dN_dxi = 0.25 * np.array([
                        [-(1-eta), -(1-xi)],
                        [(1-eta), -(1+xi)],
                        [(1+eta), (1+xi)],
                        [-(1+eta), (1-xi)]
                    ])
                    
                    # 物理坐标的导数
                    dN_dx = np.dot(dN_dxi, inv_jac)
                    
                    # 计算B^T*B
                    BTB = np.zeros((4, 4))
                    for i in range(4):
                        for j in range(4):
                            BTB[i,j] = np.dot(dN_dx[i], dN_dx[j])
                    
                    # 积分
                    Ke += w * det_jac * BTB
                
                # 单元载荷向量（假设源项为常数）
                if 'source_term' in params and pde_type == 'poisson':
                    source = params['source_term']
                    fe = 0.25 * source * np.array([1, 1, 1, 1])
                else:
                    fe = np.zeros(4)
                
            # 将单元矩阵组装到全局矩阵
                for i in range(len(elem)):
                global_i = elem[i]
                F[global_i] += fe[i]
                
                for j in range(len(elem)):
                    global_j = elem[j]
                    K[global_i, global_j] += Ke[i, j]
                F[global_i] += fe[i]
                
                for j in range(len(elem)):
                    global_j = elem[j]
                    K[global_i, global_j] += Ke[i, j]
                F[global_i] += fe[i]
                
                for j in range(len(elem)):
                    global_j = elem[j]
                    K[global_i, global_j] += Ke[i, j]
        
        # 应用边界条件
        bc_indices = []
        for i in range(num_nodes):
            x_idx = i % num_nodes_x
            y_idx = i // num_nodes_x
            
            # 检查是否在边界上
            on_boundary = False
            if x_idx == 0 and boundary_conditions['left']['type'] == 'dirichlet':
                val = boundary_conditions['left']['value']
                on_boundary = True
            elif x_idx == num_nodes_x-1 and boundary_conditions['right']['type'] == 'dirichlet':
                val = boundary_conditions['right']['value']
                on_boundary = True
            elif y_idx == 0 and boundary_conditions['bottom']['type'] == 'dirichlet':
                val = boundary_conditions['bottom']['value']
                on_boundary = True
            elif y_idx == num_nodes_y-1 and boundary_conditions['top']['type'] == 'dirichlet':
                val = boundary_conditions['top']['value']
                on_boundary = True
            
            if on_boundary:
                bc_indices.append(i)
                
                # 强制施加狄利克雷边界条件
                K[i, i] = 1e10  # 大数法
                F[i] = val * 1e10
        
        # 求解线性方程组
        U = np.linalg.solve(K, F)
        
        # 重构解矩阵
        solution = U.reshape((num_nodes_y, num_nodes_x))
        
        # 创建网格信息
        mesh = {
            "num_nodes": num_nodes,
            "num_elements": len(elements),
            "nodes": [[X[i%num_nodes_x, i//num_nodes_x], Y[i%num_nodes_x, i//num_nodes_x]] for i in range(num_nodes)],
            "elements": elements,
            "element_type": element_type
        }
        
# 此代码片段显示的是一个不完整的上下文，从错误提示来看，可能是原代码中该 `return` 语句位置有误。
# 但仅依据当前选择内容，若要让 `return` 语句能正常使用，需将其置于完整函数上下文中。
# 由于此处推测该代码是在 `finite_element` 函数中，且假定 `solution` 和 `mesh` 变量已正确定义，以下为修正后的形式

        # 以下代码是假设在 `finite_element` 函数中，补充了前面未完成的部分计算
        # 假设这里已经完成了刚度矩阵和载荷向量的组装，求解了线性方程组得到 solution
        # 同时 mesh 变量也已正确定义
        # 求解线性方程组 Kx = F
        try:
            solution = np.linalg.solve(K, F)
        except np.linalg.LinAlgError as e:
            raise ValueError(f"求解线性方程组失败: {str(e)}") from e
            # 构建网格信息
            mesh = {
                'nodes': np.column_stack([X.flatten(), Y.flatten()]).tolist(),
                'elements': elements
            }
            return {"solution": solution.tolist(), "mesh": mesh}
        except np.linalg.LinAlgError as e:
            raise ValueError(f"求解线性方程组失败: {str(e)}")

    def finite_difference_pde(self, params: Dict[str, Any]) -> Dict[str, Any]:
        """
        有限差分法求解偏微分方程
        
        参数:
            pde_type: 偏微分方程类型 ('heat' 热传导方程, 'wave' 波动方程, 'laplace' 拉普拉斯方程)
            boundary_conditions: 边界条件，格式为:
                {
                    'left': {'type': 'dirichlet', 'value': 0},
                    'right': {'type': 'dirichlet', 'value': 0},
                    'top': {'type': 'dirichlet', 'value': 0},
                    'bottom': {'type': 'dirichlet', 'value': 0}
                }
            initial_conditions: 初始条件函数（仅对抛物型和双曲型方程需要）
            x_range: x轴范围 [x_start, x_end]
            y_range: y轴范围 [y_start, y_end] (仅对椭圆型方程需要)
            t_range: 时间范围 [t_start, t_end] (仅对抛物型和双曲型方程需要)
            dx: x轴空间步长
            dy: y轴空间步长 (仅对椭圆型方程需要)
            dt: 时间步长 (仅对抛物型和双曲型方程需要)
            num_steps_x: x轴网格点数
            num_steps_y: y轴网格点数 (仅对椭圆型方程需要)
            num_steps_t: 时间步数 (仅对抛物型和双曲型方程需要)
            
        返回:
            solution: 解的二维或三维数组
        """
        pde_type = params['pde_type']
        boundary_conditions = params['boundary_conditions']
        x_range = params.get('x_range', [0, 1])
        dx = params.get('dx', 0.1)
        num_steps_x = params.get('num_steps_x', int((x_range[1] - x_range[0])/dx) + 1)
        
        if pde_type == 'heat':  # 抛物型方程（热传导）
            # 需要时间参数
            t_range = params['t_range']
            dt = params['dt']
            num_steps_t = params.get('num_steps_t', int((t_range[1] - t_range[0])/dt) + 1)
            
            # 稳定性检查
            if dt > dx**2 / 2:
                raise ValueError("Time step too large for stability")
            
            # 初始化网格
            u = np.zeros((num_steps_t, num_steps_x))
            
            # 设置初始条件
            if 'initial_conditions' in params:
                x = np.linspace(x_range[0], x_range[1], num_steps_x)
                u[0, :] = params['initial_conditions'](x)
            else:
                u[0, :] = 0
            
            # 设置边界条件
            if boundary_conditions['left']['type'] == 'dirichlet':
                u[:, 0] = boundary_conditions['left']['value']
            
            if boundary_conditions['right']['type'] == 'dirichlet':
                u[:, -1] = boundary_conditions['right']['value']
            
            # 显式有限差分法求解
            for n in range(num_steps_t - 1):
                for i in range(1, num_steps_x - 1):
                    u[n+1, i] = u[n, i] + (dt/dx**2) * (u[n, i+1] - 2*u[n, i] + u[n, i-1])
            
            return {"solution": u.tolist()}
            
        elif pde_type == 'wave':  # 双曲型方程（波动）
            # 需要时间参数和初始速度
            t_range = params['t_range']
            dt = params['dt']
            num_steps_t = params.get('num_steps_t', int((t_range[1] - t_range[0])/dt) + 1)
            
            # 稳定性检查
            if dt > dx:
                raise ValueError("Time step too large for stability")
            
            # 初始化网格
            u = np.zeros((num_steps_t, num_steps_x))
            
            # 设置初始条件
            if 'initial_conditions' in params:
                x = np.linspace(x_range[0], x_range[1], num_steps_x)
                u[0, :] = params['initial_conditions'](x)
            else:
                u[0, :] = 0
            
            # 设置初始速度
            if 'initial_velocity' in params:
                initial_velocity = params['initial_velocity']
            else:
                initial_velocity = lambda x: 0
            
            # 计算第一层
            for i in range(num_steps_x):
                u[1, i] = u[0, i] + dt * initial_velocity(x[i])
            
            # 设置边界条件
            if boundary_conditions['left']['type'] == 'dirichlet':
                u[:, 0] = boundary_conditions['left']['value']
            
            if boundary_conditions['right']['type'] == 'dirichlet':
                u[:, -1] = boundary_conditions['right']['value']
            
            # 有限差分法求解
            for n in range(1, num_steps_t - 1):
                for i in range(1, num_steps_x - 1):
                    u[n+1, i] = 2*u[n, i] - u[n-1, i] + (dt/dx)**2 * (u[n, i+1] - 2*u[n, i] + u[n, i-1])
            
            return {"solution": u.tolist()}
            
        elif pde_type == 'laplace':  # 椭圆型方程（拉普拉斯方程）
            # 需要y轴参数
            y_range = params['y_range']
            dy = params['dy']
            num_steps_y = params.get('num_steps_y', int((y_range[1] - y_range[0])/dy) + 1)
            
            # 初始化网格
            u = np.zeros((num_steps_y, num_steps_x))
            
            # 设置边界条件
            if boundary_conditions['left']['type'] == 'dirichlet':
                u[:, 0] = boundary_conditions['left']['value']
            
            if boundary_conditions['right']['type'] == 'dirichlet':
                u[:, -1] = boundary_conditions['right']['value']
            
            if boundary_conditions['top']['type'] == 'dirichlet':
                u[0, :] = boundary_conditions['top']['value']
            
            if boundary_conditions['bottom']['type'] == 'dirichlet':
                u[-1, :] = boundary_conditions['bottom']['value']
            
            # 迭代求解（雅可比迭代法）
            max_iter = 1000
            tol = 1e-4
            
            for iter in range(max_iter):
                u_old = u.copy()
                
                for i in range(1, num_steps_y - 1):
                    for j in range(1, num_steps_x - 1):
                        u[i, j] = 0.25 * (u_old[i+1, j] + u_old[i-1, j] + u_old[i, j+1] + u_old[i, j-1])
                
                # 检查收敛
                if np.max(np.abs(u - u_old)) < tol:
                    break
            
            return {"solution": u.tolist(), "iterations": iter, "converged": bool(iter < max_iter)}
            k4 = h * dy_dx(x[i] + h, y[i] + k3)
            
            y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
            k3 = h * dy_dx(x[i] + h/2, y[i] + k2/2)
